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I.     INTRODUCTION 

A.       THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a  closed  container  designed  to  transfer  a 
large  amount  of  heat  in  rotating  machinery.  Since  the  heat  pipe  operates 
on  a  closed  two-phase  cycle,  the  heat  transfer  capacity  is  greater  than  for 
solid  conductors.  Also,  the  thermal  response  time  is  less  than  with  solid 
conductors.  The  three  major  elemental  parts  of  the  rotating  heat  pipe  are: 
a  cylindrical  evaporator,  a  truncated  cone  condenser,  and  a  fixed  amount 
of  working  fluid  as  shown  in  figure  1. 

An  annulus  is  formed  by  the  working  fluid  in  the  evaporator.  This 
occurs  at  rotationary  speeds  above  the  critical  speed.  The  addition  of  heat 
to  the*  evaporator  vaporizes  the  working  fluid.  A  pressure  differential 
between  the  evaporator  and  the  condenser  causes  the  vapor  to  flow 
towards  the  condenser.  The  vapor  is  transported  to  the  condenser  with  its 
latent  heat  of  vaporization.  Condensation  of  the  vapor  on  the  inner  wall  is 
caused  by  external  cooling.  This  condensation  releases  the  latent  heat  of 
evaporation.  This  condensate  is  forced  to  flow  back  to  the  evaporator  by 
a  component,  acting  along  the  condenser  wall,  of  the  centrifugal  force 
which  is  caused  by  the  rotation  of  the  heat  pipe.  As  the  condensate 
collects  in  the  evaporator  the  cycle  is  repeated. 

Since  the  evaporator  and  condenser  portions  of  a  heat  pipe  function 
independently,  needing  only  common  liquid  and  vapor  streams,  the  area 
over  which  heat  is  introduced  can  differ  in  size  and  shape  from  the  area 
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over  which  it  is  rejected,  provided  that  the  rate  at  which  the  liquid  is 
vaporized  does  not  exceed  the  rate  at  which  it  can  be  condensed. 
Therefore,  high  heat  fluxes  generated  over  relatively  small  areas  can  be 
dissipated  over  larger  areas  with  reduced  heat  fluxes;  allowing  a 
cylindrical  evaporator  and  a  truncated  cone  condenser. 

Capillary  action  acts  to  drive  the  condensate  back  to  the  evaporator 
in  a  conventional  heat  pipe.  No  limitation  due  to  capillary  action  is 
encountered  in  a  rotating  heat  pipe  nor  are  external  pumps  or  gravity 
depended  on  for  the  flow  of  the  working  fluid.  Therefore,  the  rotating  heat 
pipe  can  be  used  in  any  orientation  [Ref.  1]. 

B.       OPERATING  LIMITS  OF  A  ROTATING  HEAT  PIPE 

The  first  theoretical  investigation  of  the  rotating  heat  pipe  conducted 
at  the  Naval  Postgraduate  School  was  performed  by  Ballback  [Ref.  2]  in 
1969.  Various  fluid  dynamic  mechanisms  limit  the  performance  of  a  rotating 
heat  pipe.  Ballback  [Ref.  2]  studied  these  mechanisms  and  an  estimation  of 
the  sonic  limit,  boiling  limit,  entrainment  limit,  and  condensing  limit  of 
performance  was  made. 

1.      The  Sonic  Limit 

The  maximum  flow  of  the  vapor  is  set  by  the  choked  flow 
condition  in  the  rotating  heat  pipe.  This  limiting  vapor  flow  rate  occurs 
when  the  heat  flux  is  increased  and  limits  the  amount  of  energy  the  vapor 
can  transport.  The  rotating  heat  pipe  effectiveness  is  reduced  due  to  this 
limitation.  The  limiting  heat  transfer  rate  due  to  this  condition  is: 


Q,  -  PVUA*A 

and  the  vapor  velocity  is  considered  to  be  sonic, 


(1) 


U-c-  (gkRD 


(2) 


where 

Uv  =  velocity  of  the  vapor  in  ft/sec,  and 

A  =  cross  sectional  area  for  the  vapor  flow  in  ft 

c  =  sonic  velocity  in  ft/ sec 

g0  =  gravitational  constant 

k  =  ratio  of  specific  heats 

R  =  gas  constant  in  ft-lbf/lbm  R 

T  =  absolute  temperature  in  degrees  Rankine 

pv  =  density  of  vapor  in  lbm/ft 

2.      The  Boiling  Limit 

The  transition  from  nucleate  to  film  boiling  was  hypothesized  by 
Kutateladze  [Ref.  3]  to  be  a  completely  hydrodynamic  process.  He 
determined  the  following  theoretical  formula  for  predicting  the  burnout 
flux: 


Qt  -  ^^^^(p^p^} 


(3) 


where 

K  =  constant  value 

Au  =  heat  transfer  area  in  the  boiler  in  ft 


hfg  =  latent  heat  vaporization  in  Btu/lbm 

o  =  surface  tension  in  lbf/ft 

g  =  acceleration  of  gravity  in  ft/hr 

p*  =  density  of  fluid  in  lbm/ft 

pv  =  density  of  vapor  in  lbm/ft  . 
A  constant  value  for  K  in  the  range  of  0.13  to  0.19  was  suggested  by  the 
experimental  data  obtained  by  Kutateladze. 

3.      The  Entrainment  Limit 

The  flooding  constraint  in  a  wickless  heat  pipe  was  examined  by 
Sakhuja  [Ref.  4].  The  correlation  he  developed  is: 


AxC\JgD(f>f-PJ)9 


(4) 


where 


Qt  =  heat  transfer  rate  in  Btu/hr 

Ax  =  flow  rate  in  ft 

C  =  dimensionless  constant,  0.725  for  tube  with  sharp  edged  flange 

hfg  =  latent  heat  of  vaporization  in  Btu/lbm 

g  =  acceleration  due  to  gravity  in  ft/hr 

D  =  inside  diameter  of  heat  pipe  in  ft 

pf  =  density  of  the  fluid  in  lbm/ft 

pv  =  density  of  the  vapor  in  lbm/ft  . 


4.        The  Condensing  Limit 

The  condenser  section  of  a  rotating  heat  pipe  was  modeled  as  a 
truncated  cone  by  Ballback  [Ref.  2].  Using  this  model,  the  condensation 
limitation  for  a  rotating  heat  pipe  was  determined  by  Ballback  [Ref.  2].  He 
developed  the  following  condensation  limit: 

Qt  .  w{2  frP/"2 W-rJVM  w^utotft-jfF  (5) 

where 

Qt.  =  total  heat  transfer  rate  in  Btu/hr 

kf  =  thermal  conductivity  of  the  condensate  film  in  Btu/hr-ft-F 

pf  =  density  of  fluid  in  lbm/ft 

o>  =  angular  velocity  in  1/hr 

hfg  =  latent  heat  of  vaporization  in  Btu/lbm 

Ts  =  saturation  temperature  in  F 

TM  =  inside  wall  temperature  in  F 

Uf  =  viscosity  of  fluid  in  lbm/ft-hr 

$  =  half  cone  angle  in  degrees 

R0  =  minimum  wall  radius  in  ft 

L  =  length  along  the  wall  of  the  condenser  in  feet. 
The    geometry    and    speed    of   the    rotating    heat    pipe,    and    the    physical 
properties  of  the  working  fluid  comprise  the  condensing  limit  equation. 

For  a  rotating  heat  pipe  with  the  physical  characteristics  as 
shown  in  Table  I,  the  amount  of  heat  that  can  be  transferred  from  the 
rotating    heat    pipe    is    limited    by    the    condensing    limit.    However,    the 


limitations  imposed  by  the  sonic  limit,  boiling  limit,  and  entrainment  limit 
may  become  important  as  the  heat  pipe  geometry  and  operating  conditions 
vary. 

TABLE  I.     ROTATING  HEAT   PIPE   SPECIFICATIONS 


Length 


9.000  inches 


Minimum  Diameter 


1.55  inches 


Wall  Thickness 


0.03125  inches 


Internal  Half  Angle 


1.000  degree 


Rotating  Speed 


3600  RPM 


To  enhance  the  heat  transfer  capacity  of  the  rotating  heat  pipe, 
internally  finned  condensers  have  been  used  to  raise  the  condensing  limit 
line.  Thinner  films  occur  near  the  ridges  of  the  fins  while  thicker  films 
occur  in  the  troughs.  The  thinner  film  on  the  ridges  provides  a  lower 
thermal  resistance  to  heat  flow,  while  the  thicker  film  in  the  trough 
provides  a  higher  resistance.  A  compromise  between  the  improvement  on 
the  ridges  and  the  degradation  in  the  troughs  is  necessary  for  an  overall 
heat  transfer  improvement  [Ref.  5]. 


C.       ANALYSIS  OP  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 

Schafer  [Ref.  6]  developed  an  analytical  model  for  a  heat  pipe  with 
a  triangular  fin  profile  (figure  2).  This  model  was  developed  in  order  to 
raise  the  condensing  limit  by  the  addition  of  internal  fins.  An  assumption 
of  one-dimensional  heat  conduction  through  the  wall  and  fin  was  made  for 
Schafer' s  model. 

A  two-dimensional  heat  conduction  model  using  a  Finite  Element  Method 
was    developed     for    this     same    case    by     Corley     [Ref.     7].     A     parabolic 


temperature  distribution  along  the  fin  surface  was  assumed  by  Corley  [Ref. 
7].  A  significant  improvement  in  the  predicted  heat  transfer  performance 
was  indicated  by  his  results.  By  using  the  two-dimensional  model  an 
increase  of  approximately  75  percent  in  the  heat  transfer  performance  was 
seen  over  the  results  from  the  use  of  the  one  dimensional  model.  However, 
Corley  [Ref.  7]  noted  that  at  the  fin  apex  an  error  as  great  as  50  percent 
was  possible  which  could  result  in  the  total  heat  transfer  being  in  error 
as  much  as  15  percent. 

A  modification  was  made  to  Corley's  computer  program  by  Tantrakul 
[Ref.  8].  In  order  to  minimize  the  heat  transfer  error  at  the  apex  of  the 
fin  Tantrakul  increased  the  number  of  finite  elements  used.  His  results 
with  this  modification  converged  with  the  results  of  Corley. 

Purnomo  developed  a  linear  triangular  finite  element  model  (figure  3) 
used  in  a  two-dimensional  Finite  Element  Method  solution.  Purnomo's  [Ref. 
1]  Finite  Element  Method  program  also  worked  and  converged.  To  maximize 
the  heat  transfer  from  the  rotating  heat  pipe  the  condenser  geometry  was 
varied.  Using  Purnomo's  code  parametric  studies  were  conducted.  However, 
the  best  geometry  was  not  indicated  in  these  studies.  Purnomo's  code  was 
written  to  perform  one  analysis  at  a  time.  Davis  [Ref.  9]  modified 
Purnomo's  code  to  allow  for  numerous  analysis  to  be  made  using  the 
optimization  code  COPES/CONMIN.  Davis'  Finite  Element  .  Method  code 
incorporating  the  optimization  worked  and  converged,  resulting  in  an 
optimum  design  for  an  internally  finned  rotating  heat  pipe. 


Figure  2.        Internally    Finned     Condenser    Geometry,     Showing     Fins, 
Troughs  and  Lines  of  Symmetry. 


e  •  thicKnasa 

b  •  heightt  of  th«  fin 

9     ■  condenser  half  angle 
a   »  half  of  th«  fin  width 
c»  half  of  tarn   trough  vidth 

Figure  3.     Condenser  Geometry   Considered  with  40  Linear  Triangular 
Finite  Elements. 


D.        THESIS  OBJECTIVES 

The  objectives  of  this  thesis  were: 

1.  To  modify  Davis'  [Ref.  9]  computer  program  so  that  it  is  compatible 
with  the  ADS  (Automated  Design  Synthesis)  program  [Ref.  10]  and 
can  be  used  for  analysis  and  automated  design  of  rotating  heat 
pipes. 
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2.  To  use  the  resulting  program  to  obtain  an  optimum  design  for  an 
internally  finned  rotating  heat  pipe  to  obtain  experimental  data  to 
compare  with  the  analytical  results. 

3.  To  use  the  resulting  program  to  obtain  numerical  results  in  place  of 
data  obtained  from  costly  experimental  operations. 
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II.     NUMERICAL  OPTIMIZATION 

A.        BACKGROUND 

The  parameter  that  is  minimized  or  maximized  during  the  design 
process  is  called  the  design  objective.  The  design  objective  is  minimized 
or  maximized  by  changing  the  design  variables  within  the  design  constraint 
limitations.  This  process  is  called  numerical  optimization.  An  assortment  of 
physical,  aesthetic,  economic  and,  on  occasion,  political  limitations  must  be 
met  by  the  design  constraints  for  the  design  to  be  acceptable.  For  the 
optimization  process  to  work,  the  design  criteria  must  be  described  in 
numerical  terms.  This  is  not  always  easy. 

A  computer  program  can  be  written  to  perform  tedious  and  repetitive 
calculations  necessary  to  optimize  the  problem  once  it  is  stated  in 
numerical  terms.  For  this  reason,  computer  analysis  is  commonplace  in  most 
engineering  organizations.  For  example,  in  heat  transfer  design  the 
configuration,  materials,  and  method  of  heat  removal  may  be  defined  and 
a  finite  element  analysis  computer  code  is  used  to  calculate  temperatures, 
heat  transfer  rates,  and  other  response  quantities  of  interest.  If  any  of 
these  parameters  are  not  within  prescribed  bounds,  the  engineer  may 
change  the  method  of  cooling  or  other  defined  quantity  and  rerun  the 
program.  The  engineer  makes  the  actual  design  decisions,  the  computer 
code  only  provides  the  analysis  of  a  proposed  design.  This  is  the  commonly 
used  approach  which  is  called  computer-aided  design. 
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Analysis  codes  are  commonly  used  for  tradeoff  studies.  For  example, 
an  analysis  code  might  be  run  on  the  distance  a  truck  can  go  on  a  tank 
of  fuel.  For  different  loads,  different  distances  are  calculated  which  can 
be  used  in  a  range-payload  study. 

Fully  automated  design  is  the  logical  next  step  to  computer-aided 
design.  The  computer  makes  the  actual  design  decisions  or  trade-off 
studies  based  on  input  criteria  in  fully  automated  design.  Minimal 
information  is  requested  from  the  operator  during  the  actual  design 
process.  Numerical  optimization  offers  numerous  improvements  over  the 
traditional  approach  to  design.  These  improvements  include:  time  reduction 
in  design  decision  making;  a  rational,  directed  design  procedure;  and  the 
procedure  is  unbiased  by  intuition  or  experience.  The  probability  of 
obtaining  a  non-traditional  solution  is  thereby  improved.  Engineering 
intuition  and  experience  are  still  necessary  to  decide  if  the  design  obtained 
is  an  improvement  and  feasible. 

B.       AUTOMATED  DESIGN   SYNTHESIS  (ADS) 

Vanderplaats  [Ref.  10]  developed  a  general  purpose  numerical 
optimization  program  containing  a  variety  of  algorithms,  ADS.  ADS  is  a 
FORTRAN  program  that  optimizes  a  numerically  defined  objective  function 
subject  to  a  set  of  constraint  limits.  The  solution  of  the  problem  is 
separated  into  three  levels: 

1.  Strategy    -    Optimization    strategy     such    as    Augmented     Lagrange 
Multiplier  method  or  Sequential  Linear  Programming. 

2.  Optimizer  -  Actual  algorithm  to  perform  the  optimization . 

3.  One-Dimensional  Search  -  Line  search  routine  used  by  optimizer. 
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Flexibility  to  solve  a  wide  variety  of  engineering  design  problems  is 
given  by  the  combinations  of  nine  strategies,  five  optimizers,  and  eight 
one-dimensional  search  options.  The  following  definitions  are  necessary  to 
discuss  the  use  of  ADS: 


1.  Design  Variables  -  Those  parameters  which  the  optimization  program 
is  permitted  to  change  within  allowed  bounds  in  order  to  improve 
the  design.  Design  variables  appear  only  on  the  right  hand  side  of 
an  equation  and  are  continuous. 

2.  Design  Constraints  -  An  inequality  constraint  requires  that  some 
function  of  the  design  variable(s)  remain  less  than  a  specified  value. 
Design  constraints  may  be  linear  or  nonlinear,  implicit  or  explicit, 
but  they  must  be  continuous  functions  of  the  design  variable. 

3.  Objective  Function  -  The  parameter  which  is  going  to  be  minimized 
or  maximized  during  the  optimization  process.  The  objective  function 
may  be  linear  or  nonlinear,  implicit  or  explicit,  and  must  be  a 
continuous  function  of  the  design  variables.  The  objective  function 
usually  appear  on  the  left  side  of  an  equation. 


C.        PROGRAMMING  GUIDELINES 

Any  computer  code  developed  for  engineering  analysis  should  be 
written  in  such  a  way  that  it  is  easily  coupled  to  a  general  purpose 
optimization  program  such  as  ADS.  Therefore,  a  general  programming 
practice  is  outlined  here  which  in  no  way  inhibits  the  use  of  the  computer 
program  in  its  traditional  role  as  an  analytical  tool,  but  allows  for  simple 
adaption  to  ADS. 

ADS  is  called  by  a  user-supplied  calling  program.  ADS  does  not  call 
any  user-supplied  subroutines.  Instead,  ADS  returns  control  to  the  calling 
program  when  function  or  gradient  information  is  needed.  The  required 
information  is  evaluated  and  ADS  is  called  again.  This  provides  considerable 
flexibility  in  program  organization  and  restart  capabilities.  Various  internal 
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parameters  are  defined  on  the  first  call  to  ADS  which  work  well  for  the 
"average"  optimization  task.  However,  it  is  often  desirable  to  change  these 
in  order  to  gain  maximum  utility  of  the  program.  Figure  4  is  the  program 
flow  diagram  for  the  case  where  the  user  wishes  to  over-ride  one  or  more 
internal  parameters,  such  as  scaling,  convergence  criteria,  or  maximum 
number  of  iterations. 

After  initialization  of  basic  parameters  and  arrays,  the  information 
parameter,  INFO,  is  set  to  -2.  ADS  is  then  called  to  initialize  all  internal 
parameters  and  allocate  storage  space  for  internal  arrays.  Control  is  then 
returned  to  the  user,  at  which  point  these  parameters,  for  example 
convergence  criteria,  can  be  overridden  if  desired.  At  this  point,  INFO  will 
have  a  value  of  1  and  the  user  must  evaluate  the  objective  function,  OBJ, 
and  constraint  functions.  ADS  is  called  again  and  the  optimization  proceeds. 
Since,  in  this  case,  the  gradient  calculation  control,  IGRAD,  has  a  value  of 
zero,  all  gradient  information  is  calculated  by  finite  difference  within  ADS. 
When  INFO  has  a  value  of  zero,  optimization  is  complete. 
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BEGIN 

DIMENSION   ARRAYS 

DEFINE  BASIC  VARIABLES 

IGRAD=0  (USE   FINITE  DIFFERENCE  GRADIENTS) 

INFO  =   -2 

CALL  ADS  (INFO...) 

IF  INFO  =  0,  EXIT.  ERROR  WAS   DETECTED 

ELSE 

OVER-RIDE  DEFAULT   PARAMETERS  IN  ARRAYS  WK  AND  IWK  IF 

DESIRED 

CALL  ADS   (INFO...) 

NO  YES 

INFO  =  0 


EVALUATE  OBJECTIVE  EXIT  OPTIMIZATION 

AND  CONSTRAINTS  IS   COMPLETE 


Figure    4.        Program    Flow    Logic:    Over-Ride    Default    Parameters,    Finite 
Difference  Gradients  [Ref.  10]. 
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m.     FINITE  ELEMENT'   SOLUTION 

A.        REVIEW   OF  THE   PREVIOUS  ANALYSIS 

As  stated  previously,  the  heat  transfer  solution  for  a  one-dimensional 
model  of  an  internally  finned  rotating  heat  pipe  was  studied  by  Schafer 
[Ref.  6].  The  two-dimensional  model  was  studied  by  Corley  [Ref.  7].  The 
same  assumptions  and  boundary  conditions,  similar  to  those  used  in  the 
Nusslet  analysis  of  film  condensation  on  a  flat  wall,  and  based  upon  the 
analysis  of  Ballback  [Ref.  2]  were  used  for  both.  The  more  important  of 
these  assumptions  are: 

1.  steady  state  operation, 

2.  film  condensation,  as  opposed  to  drop  wise  condensation, 

3.  laminar    flow    of   the    condensate    film    along    both    the    fin    and    the 
trough, 

4.  static  balance  of  forces  within  the  condensate, 

5.  one-dimensional  conduction  heat  transfer  through  the  film  thickness 
(no  convective  heat  transfer  in  the  condensate  film), 

6.  no  liquid-vapor  interfacial  shear  forces, 

7.  no  condensate  subcooling, 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of  the  wall  section 
(symmetry  conditions),  as  shown  in  figure  5, 

9.  saturation  temperature  at  the  fin  apex, 

10.  zero  film  thickness  at  the  fin  apex,  and 

11.  negligible  curvature  of  the  condenser  wall. 
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Figure  3  shows  the  linear  triangular  finite  element  model  developed 
by  Purnomo  [Ref .  1]  for  use  in  obtaining  a  two-dimensional  Finite  Element 
solution. 

The  assumption  that  was  used  by  Corley  [Ref.  7]  that  the  fin  apex  was 
at  the  saturation  temperature  of  the  working  fluid  was  modified  by 
Purnomo  [Ref.  1].  The  value  of  the  temperature  at  the  apex  of  the  fin  was 
allowed  to  float  and  a  parabolic  temperature  distribution  was  assumed  along 
the  fin  surface. 

Purnomo's  problem  statement  for  the  formulation  of  the  Finite  Element 
Method  as  shown  in  figure  5  is: 


dx2  dy- 

with  the  following  -boundary  conditions: 

1.  along  boundary  1,      -K  dT/dn-h^T-Tsat) 

2.  along  boundary  3,      -K  dTldn-h^T-T™) 

dT  n 

3.  along  boundaries  2  and  4,      -0 

dn 

A   detailed   description    of  the    numerical  formulation  is    presented   in    his 
thesis. 

Davis     used     Constrained     Function     Minimization     (CONMIN)     as     an 
optimization  program.  CONMIN   is  a   FORTRAN  program   in  subroutine  form. 
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b.c. 


3*T  *  32T 
3x2    3y2 


3T 

a)  -  :<  g-jj  =  hL  CT  -  Tgat)  Along  Boundary  Cl] 

b)  -  k  t^  =  h2  (T  -  T^)  Along  Boundary  [3] 

3  T 

c)  =rr  =  3  Along  3oundarias  C2]  and  £«*.] 


FLq^e  «5"   ,  Differential  Equation  and  Boundary   Conditions   Considered 
ui  the  Analysis  of  Purnomo  [Ref.  i]. 
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Vanderplaats  [Ref.  11]  developed  the  Control  Program  for  Engineering 
Synthesis  (COPES)  as  a  main  program  to  simplify  the  use  of  CONMIN.  Davis' 
computer  program  was  written  in  subroutine  form  with  SUBROUTINE  ANALIZ 
(ICALC)  as  the  main  routine.  The  name  ANALIZ  is  compatible  with  the  COPES 
program  and  ICALC  is  a  calculation  control.  Subroutine  ANALIZ  calls  other 
subroutines  as  needed: 


1.  the  routine  "COORD"   used  to   define  positions  of  system   coordinate 
points 

2.  the  routine  "FORMAF"   used  to  formulate  the  Finite  Element  Method 
equations, 

3.  the  routine  "BANDEC"  as  an  equation  solver  for  a  symmetric  matrix 
which  has  been  transformed  into  banded  form,  and 

4.  the  routine  "DPLORT"  used  to  compute  the  roots  of  a  real  polynomial 
using  a  Newton-Rap hson  derivative  technique. 


B.       THE  COMPUTER  PROGRAM  DEVELOPMENT 

The  basis  for  the  present  analysis  code  is  Davis'  [Ref.  9]  two- 
dimensional  finite  element  program.  Davis'  code  was  checked  for  validity  as 
the  first  task  undertaken  in  the  development  of  this  thesis.  An  error  was 
discovered  in  calculating  the  fin  condensate  film  thickness  (AZS).  In  the 
initial  calculation  of  HDEN  only,  the  cubed  term  was  merely  multiplied  by 
three.  The  correct  form  of  the  equation  is  shown  below: 

HDEN  -  -o^P-b^P+ziT^-TJ 

The  effect   of   this    error    was   minimal    since    subsequent   equations    were 
correct,  a  0.00016%  difference  in  the  condensate  level  was  noted. 
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The  next  task  undertaken  was  to  adapt  the  analysis  code  to  permit 
automated  design  analysis  using  ADS.  Many  modifications  were  made,  some 
of  which  are  mentioned  here.  The  program  was  rewritten  to  include  a  main 
program  from  which  ADS  is  called  and  subroutines  to  perform  various 
mathematical  functions.  Subroutine  ANALIZ  was  deleted  since  the  double 
precision  version  of  ADS  (DADS)  was  used.  Modifications  were  also  made  to 
generalize  the  code  and  minimize  the  changes  needed  when  varying  the 
number  of  finite  elements  used. 

A  listing  of  the  revised  computer  program  is  included  as  the 
Appendix. 

C.       DESIGN  OPTIMIZATION 

There  are  thirteen  parameters  that  can  be  used  as  design  variables. 
There  are  geometric  or  functional  parameters  of  the  rotating  heat  pipe  or 
the  properties  of  the  working  fluid  or  environment.  The  possible  design 
variables,  possible  constraint  functions,  and  the  objective  function  are 
listed  below  in  Fortran.  This  code  can  pursue  a  wide  variety  of  design 
problems. 

The  addition  of  fins  by  the  designer  increases  the  surface  area  which 
increases  the  heat  transfer  rate  through  the  condenser  wall.  However,  the 
addition  of  fins  decreases  the  cross-sectional  area  through  each  fin  for 
conduction  and  decreases  the  trough  width  which  increases  the  condensate 
thickness  in  the  trough.  The  increased  condensate  thickness  decreases  the 
heat  transfer  and  if  increased  to  the  point  of  covering  the  fins  it  could 
dramatically  reduce  the  heat  transfer  through  the  fin. 
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TABLE  H.  DESIGN  VARIABLES 


DESIGN  VARIABLES 


BFIN  (fin  height) 

CANGL   (cone  half  angle) 

CLI  (condenser  length 

FANGL  (fin  half  angle) 

HINF  (external  convective  heat  transfer  coefficient) 

NFIN  (number  of  fins) 

R2I  (intermediate  radius) 

RBASEI  (inside  radius  of  condenser  base) 

RMP  (rotational  speed  of  the  heat  pipe) 

THICKI  (condenser  wall  thickness) 

TINF  (external  temperature) 

TSS  (saturation  temperature) 

TZ  (nodal  point  temperature) 


CONSTRAINT  FUNCTIONS 


BOA  (ratio  of  fin  height  over  fin  width) 
ZOA  (ratio  of  trough  width  over  fin  width) 
DEL(NI)  (condensate  thickness) 


OBJECTIVE  FUNCTION 


OBJ 


The  purpose  of  the  design  study  was  to  determine  the  fin  height, 
number  of  fins,  and  fin  half  angle  which  would  maximize  the  heat  transfer 
rate.  It  was  then  decided  that  the  design  variables  would  be  BFIN,  FANGL, 
and  NFIN.  The  number  of  fins  was  chosen  vice  the  ratio  of  trough  width 
to  fin  width  (ZOA)"  as  was  previously  used.  This  decision  was  made  since 
it  is  easier  to'  think  in  terms  of  the  number  of  fins  vice  a  ratio.  Other 
potential  design  variables  were  held  constant.  The  objective  function  to  be 
maximized  was  OBJ=QT+QTF,  the  heat  transfer  through  the  fin  plus  the  heat 
transfer  through  the  trough. 

The  code  was  run  with  the  three  design  variables  using  the  internal 
scalar  default  parameters  in  ADS.  The  objective  function  was  calculated 
using   the  input   values   of   the   design    variables.    The   ADS    program    then 
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changed  the  first  design  variable  keeping  the  other  two  design  variables 
at  the  input  values.  The  calculations  were  made  aging  yielding  a  different 
value  for  the  objective  function.  The  new  value  would  then  be  compared 
to  the  previous  value  of  the  objective  function.  If  the  difference  in  the 
objective  function  was  not  greater  then  the  internal  scalar  default  value, 
the  first  design  variable  would  be  returned  to  its  original  value  and  the 
process  repeated  with  the  second  design  variable.  If  the  difference  in  the 
objective  function  was  still  not  crreater  then  the  internal  scalar  default 
value,  the  second  design  variable  would  be  returned  to  its  initial  value  and 
the  process  repeated  with  the  third  design  variable. 

Each  time  the  program  was  run,  the  optimization  code  would  choose 
BFIN  as  the  design  variable  to  change  first  as  it  had  the  greater  influence 
on  the  objective  function.  The  remaining  two  variables  would  then  either 
be  kept  constant  or  the  number  of  fins  would  be  maximized  and  the  fin 
half  angle  minimized.  Consistent  results  were  not  obtained  with  this  method. 

To  improve  the  results,  the  internal  scalar  parameters  were  modified. 
These  modifications  included  the  constraint  tolerances,  the  absolute  and 
relative  convergence  criteria,  the  absolute  and  relative  change  in  the 
design  variables,  the  absolute  and  relative  change  in  the  objective 
function,  the  minimum  absolute  value  of  the  finite  difference  step  when 
calculation  gradients,  and  the  initial  relative  move  limit.  Better  results  were 
obtained  as  seen  in  the  higher  value  for  the  objective  function.  However, 
the  results  were  still  not  consistent  depending  on  the  initial  values  used 
for  the  three  design  variables.  At  this  point,  it  was  decided  to  concentrate 
on  one  design  variable  and  on  the  basis  of  the  previous  calculations  the 
design  variable  chosen   was  BFIN. 
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The  external  surface  temperature  was  set  equal  to  the  working  fluid 
saturation  temperature  and  the  theoretical  upper  limit  on  the  heat  transfer 
was  calculated  for  comparison.  An  assumption  is  made  that  there  is  no 
thermal  resistance  across  the  condensate  or  the  condenser  wall.  The  upper 
limit  of  the  heat  transfer  rate  was  predicted  to  be  69,492  BTU/HR  using  the 
following  formula: 


Q^  -  2Wrt(T.-rj  cv) 


•max  v    wall 

where 

h  =  outside  convective  heat  transfer  coefficient  (5000  BTU/HR '  FT    F) 

r  =  average  outside  radius  of  condenser  wall  (0.07373  FT) 

1  =  condenser  length  (0.75  FT) 

Twall  =  temperature  of  the  outside  wall  (100 *F) 

* 

Tffl  =  ambient  temperature  (60 *F) 
Based  on  engineering  judgement  certain  constraints  were  placed  on  the 
design.  These  constraints  were  applied  to  the  number  of  fins  (not  to 
exceed  400)  and  the  minimum  fin  half  angle  (not  to  be  less  than  10 
degrees).  The  values  were  based  on  structural  and  manufacturing 
considerations. 
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IV.     RESULTS 

A.  INTRODUCTION 

The  purpose  of  the  design  optimization  was  to  maximize  the  heat 
transfer  rate.  This  was  accomplished  by  using  the  computer  code  in 
conjunction  with  ADS.  Numerical  results  are  discussed  below. 

B.  CONSTRAINED  OPTIMIZATION 

In  the  design  problem  undertaken  to  determine  the  optimum  internal 
geometry  for  the  maximum  heat  transfer,  numerous  runs  were  made  for  a 
condenser  made  of  copper.  This  material  has  a  thermal  conductivity  of  230 
BTU/HR'FTT.  The  working  fluid  was  water. 

Since  the  fin  height  (BFIN)  was  the  design  variable,  the  initial  runs 
investigated  whether  there  was  an  optimum  fin  height.  Initially  the  fin  half 
angle  was  held  constant  at  10  degrees  and  the  number  of  fins  was  varied 
from  150  to  400.  In  each  case,  for  the  optimum  design,  the  fin  height  was 
maximized  and  the  trough  was  eliminated,  as  seen  in  figures  6  and  7. 

The  number  of  fins  was  then  held  constant  and  the  fin  half  angle  was 
varied  from  10  to  25  degrees  with  the  fin  height  remaining  the  design 
variable.  Once  again,  the  greatest  heat  transfer  rate  was  achieved  with  the 
highest  fin  height  for  each  number  of  fins  (figure  8.) 

As  seen  in  figure  9,  the  highest  heat  transfer  rate  achieved  was  for 
400  fins  with  a  10  degree  fin  half  angle  and  a  fin  height,  of  0.0345  inches. 
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Figure  6.     Heat  Transfer  Rate  vs.  Fin  Height. 
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Figure  7.     Heat  Transfer  Rate  vs.  Fin  Height. 
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Figure  8.    Heat  Transfer  Rate  vs.  Fin  Half  Angle  (Optimum  Fin  Height). 
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Figure  9.     Heat  Transfer  Rate  vs.  Number  of  Fins  (Optimum  Design). 
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Figure  10  shows  a  plot  of  heat  transfer  rate  versus  fin  half  angle  for 
a  condenser  with  between  100  and  400  fins.  The  heat  transfer  rate,  as  a 
function  of  fin  half  angle  for  a  constant  fin  height,  increases  as  the  fin 
half  angle  increases.  Davis  [Ref.  9]  concluded  that  the  heat  transfer  rate 
increased  with  an  increase  in  fin  half  angle.  This  is  correct  if  the  fin 
height  is  kept  constant.  As  the  fin  half  angle  increases,  the  surface  area 
of  the  fin  increases.  The  added  surface  area  also  has  a  thinner  film  of 
condensate  on  it  which  offers  lower  thermal  resistance.  The  trough  area 
decreases  and  the  condensate  film  in  the  trough  thickens,  increasing  the 
thermal  resistance.  This  degradation  does  not  offset  the  gain  in  the  heat 
transfer  rate  caused  by  the  fin. 

For  external  heat  transfer  coefficients  varying  from  1000-50,000 
BTU/HR«FT  *F,  the  same  optimum  design  geometry  for  a  maximum  heat 
transfer  rate  was  obtained,  which  is  stated  in  table  III  below.  Figure  11 
shows  the  strong  influence  the  external  heat  transfer  coefficient  has  on 
the  heat  transfer  rate. 

In  figure  12,  the  effect  of  the  rotating  speed  on  the  heat  transfer 
rate  is  seen.  As  the  rotational  speed  increases,  the  heat  transfer  rate 
increases,  this  is  caused  by  an  increase  in  the  element  heat  transfer 
coefficient  and  a  decrease  in  the  condensate  thickness  on  the  fin  which 
lowers  the  thermal  resistance. 
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Figure  11.     Heat  Transfer  Rate  vs.  Heat  Transfer  Coefficient. 
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Figure  12.     Heat  Transfer  Rate  vs.  Number  of  Fins  (RPM  Variation). 
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Figure  13.     Condensate  Level  vs.  Position  (100-400  Fins). 


34 


When  the  number  of  fins  was  increased  from  100  to  400,  maintaining 
the  same  fin  half  angle,  the  condensate  level  decreased  with  the  increased 
heat  transfer  rate  (figure  13).  This  occurred  because  the  thinner  film  over 
the  fins  decreased  the  resistance  across  the  film  which  raised  the 
temperatures  along  the  fm,  which  in  combination  with  the  lower  height  fins 
increased  the  temperature  on  the  outside  of  the  pipe.  This  increase  in 
temperature  brings  the  operation  closer  to  the  condensing  limit.  When  the 
fin  half  angle  was  increased  for  a  specified  number  and  height  of  f\ns,  the 
condensate  level  decreased  due  to  the  increased  trough  width  (figure  14). 


TABLE  m.     OPTIMIZATION  RESULTS  FOR 
VARIOUS  HEAT  TRANSFER  COEFFICIENTS 


OPTIMUM  DESIGN 


Fin  Height 
Fin  Half  Angle 
Number  of  Fins 


0.0345  inches 

10.0  degrees 

400 


"external 
(BTU/HR«FT2«F) 


HEAT  TRANSFER  RATE 


(BTU/HR) 


1000 

5000 

50000 


13765 
62252 
261003 


Figures  15  and  16  show  the  effect  the  increase  of  surface  area  has 
on  the  heat  transfer  rate.  In  figure  15,  the  fin  height  for  400  fins  with  a 
10  degree  half  angle  is  plotted  against  the  heat  transfer  rate.  As  the  fin 
height  is  increased  up  to  a  maximum  value  of  0.0345  inches  the  resulting 
design  is  a  sawtooth.  Figure  16  shows  the  effect  adding  more  fins  has  on 
the  heat  transfer  rate  for  a  constant  height  fin.  The  greatest  increase  in 
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the  heat  transfer  rate  is  seen  from  the  addition  of  100  fins  from  a  smooth 
tube. 

In  figure  17,  the  ratio  of  the  actual  surface  area  over  the  surface 
area  for  a  smooth  tube  is  plotted  versus  the  number  of  fins.  The  fin 
height  and  the  fin  half  angle  are  both  held  constant.  As  expected,  the 
surface  area  ratio  increases  in  a  linear  manner  as  the  number  of  fins 
increases.  The  ratio  of  the  heat  transfer  rate  over  the  heat  transfer  rate 
for  a  smooth  tube  is  plotted  for  two  different  heat  tr?nsfer  coefficients, 
h=1000  BTU/HR*FT2»F  and  5000  BTU/HR*FT2«F.  An  increase  in  the  number 
of  fins  results  in  not  only  an  increase  in  the  area  ratio  but  also  an 
increase  in  the  heat  transfer  ratio.  The  increase  in  the  heat  transfer  ratio 
is  greatest  when  going  from  a  smooth  tube  to  a  tube  with  100  fins.  The 
heat  transfer  ratio  increase  is  greater  for  the  heat  transfer  coefficient 
equal  to  5000  BTU/HR*FT  »P.  This  is  because  the  heat  transfer  coefficient 
has  a  direct  effect  on  the  heat  transfer  rate,  that  is, 


Q  -  W-r^ 


Both  curves  approach  an  asymptotic  value.  However,  the  curve  with  the 
lower  heat  transfer  coefficient  approaches  this  asymptotic  value  with  a 
fewer  number  of  fins.  Additionally,  in  view  of  the  relatively  small  increase 
in  the  heat  transfer  ratio  by  the  addition  of  fins  for  the  lower  heat 
transfer  coefficient,  consideration  should  be  given  to  the  cost  of 
manufacturing  the  fins  versus  the  benefit  derived  by  their  addition. 
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2  3*36 

FOSiTION  (INCHES) 


LEGEND 
a  s  FIN  ANGLE  =  10 
o  =  FIN  ANGLE  =  15 
a  =  FIN  ANGLE  =  20 
+  =  FiN  ANGLE  =  25 


Figure  14.    Condensate  Level  vs.  Position  (400  Fins,  10-25  Degree  Fin 
Half  Angle). 
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Figure  15.     Heat  Transfer  Rate  vs.  Pin  Height  (400  Pins), 


38 


OS 


z 
< 

QC 


< 
X 


65000 


60000- 


55000  - 


50000" 


3 

I 

OS 


QC     45000" 


40000- 


35000- 


30000- 


25000 

0 


50      100     150     200     2S0 

NUMBER  OF  FINS 


300  350  400 


Figure  16.      Heat   Transfer   Rate    vs.    Number   of   Fins    (Constant    Fin 
Height). 
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50  100  150  200  250  300  350  400 

NUMBER  OF  FINS 


LEGEND 
a  =  H=5000 
o  =  H=1000 


Figure  17.     Heat  Transfer  and  Area  Ratios  vs.  Number  of  Fins. 


40 


C.       AUTOMATED  DESIGN   SYNTHESIS  (ADS) 

This  optimization  project  was  done  on  the  IBM  mainframe  using  the 
ADS  optimization  program.  As  stated  previously,  ADS  is  a  general  purpose 
numerical  optimization  program  with  a  variety  of  algorithms  that  can  be 
used  to  tailor  the  solution.  The  solution  is  separated  into  three  levels  in 
ADS:  strategy,  optimizer,  and  one-dimensional  search.  For  this  problem  the 
following  combination  of  algorithms  were  used: 

1.  Strategy:  Sequential  Linear  Programming 

2.  Optimizer:   Modified    Method   of   Feasible    Directions   for   constrained 
minimization 

3.  One-Dimensional     Search:     Golden      Section      Method     followed      by 
polynomial  interpolation 

The  strategy  used  linearizes  a  nonlinear  problem  by  a  first  order 
Taylor  series  expansion  of  the  objective  and  constraint  functions.  The 
solution  to  this  linear  approximation  is  obtained.  The  problem  is  linearized 
again  about  this  point  and  the  new  problem  is  solved  with  the  process 
being  repeated  until  a  precise  solution  is  achieved. 

The  optimizer  chosen  is  used  to  find  a  search  direction  which  will 
minimize  the  objective  function  while  maintaining  feasibility. 

The  combination  of  strategy,  optimizer,  and  one-dimensional  search 
chosen  is  not  the  only  one  available,  nor  is  it  necessarily  the  most 
efficient.  It  did  yield  results  that  were  maximized  and  were  within  the 
constraint  tolerances. 

The  ADS  optimization  program  is  complicated  by  the  numerous  internal 
parameters     which     must     be     changed     to     obtain     an     optimal     design. 
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Complications  arise  when  there  is  a  vast  difference  in  the  scales  of  the 
design  variables.  The  design  variables  themselves  must  be  scaled  which  is 
further  complicated  when  the  variable  itself  covers  a  wide  range.  ADS  also, 
in  certain  cases,  allows  for  constraints  to  be  violated.  In  some  instances 
this  might  be  acceptable  but  not  in  this  case. 

ADS  does  not  have  a  scoping  mechanism,  that  is  the  ability  to 
decrease  the  rate  of  change  of  the  design  variable,  and  therefore  to  obtain 
a  precise  answer  the  internal  parameters  must  be  changed  repeatedly.  ADS 
also  does  not  recognize  integers,  all  numbers  are  real  therefore,  depending 
on  the  answer  given,  it  may  be  necessary  to  round  up  or  down. 


42 


V.     CONCLUSIONS 

1.  For  an  independent  increase  in  fin  half  angle,  rotational  speed  or 
number  of  fins  an  increase  is  seen  in  the  heat  transfer  rate.  As  the 
parameters  increase,  the  heat  transfer  rate  levels  off  at  the  theoretical 
maximum  heat  transfer  rate  for  the  heat  pipe.  A  decrease  in  the  fin  half 
angle  with  a  corresponding  increase  in  the  fin  height  increases  the  heat 
transfer  rate.  If  the  fin  half  angle  is  decreased  while  the  fin  height  is 
kept  constant,  then  the  heat  transfer  rate  decreases. 

2.  Maximum  heat  transfer  occurs  for  the  same  fin  geometry 
regardless  of  the  external  heat  transfer  coefficient.  For  a  specific 
condenser  radius,  as  many  fins  as  possible  should  be  machined  with  a 
minimum  fin  half  angle  at  the  maximum  fin  height. 

3.  The  computer  code  can  be  used  for  single  analysis  or  the 
automated  design  of  an  internally  finned  rotating  heat  pipe. 

4.  The  benefit  of  adding  fins  is  dependent  on  the  external  heat 
transfer  coefficient.  Consideration  of  the  cost  of  manufacturing  the  fins 
versus  the  increase  in  the  heat  transfer  rate  should  be  made. 
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VI.     RECOMMENDATIONS 

1.  Analyze  different  shaped  fins  including  rectangular  and  curved. 

2.  Modify  the  code  to  allow  for  silmutaneous  variations  of  more  than 
one  variable. 

3.  Use  different  working  fluids  and  heat  pipe  materials  to  see  if  a 
different  internal  geometry  occurs  for  the  maximum  heat  transfer  rate. 

4.  Modify  the  code  to  use  the  DOT  optimization  program  vice  ADS. 
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APPENDIX:     PROGRAM  LISTING 


PROGRAM    OLSON 
******************************************************** 

*  * 

*  ANALYSIS  OF  ROTATING  HEAT  PIPE  ,  USING  TRIANGULAR  * 

*  ELEMENT  MODEL  * 

*  COMPILED  BY  MAJOR  IGNATIUS . S . PURNOMO  IN  JUNE  1978  * 

*  * 

*  MODIFIED  TO  PERMIT  NUMERICAL  OPTIMIZATION  * 

*  USING  COPES/CONMIN  * 

*  BY  LCDR  WILLIAM  A.  DAVIS,  JR.  IN  SEPTEMBER  1980  * 

*  MODIFIED  TO  PERMIT  USE  OF  THE  ADS  OPTIMIZATION  * 

*  ROUTINE  BY  LT.  G.L.  OLSON   IN  JUNE  1992  * 

* 

*  * 
******************************************************** 


CHARACTER* 20  NAME 

COMMON/ADS/DF(21) ,G(10) ,IDG(100) , IGRAD, INFO, IOPT, IONED, IPRINT, 
:ISTRAT,IWK(2000) ,IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W(21,30),WK(5000) 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B(3) , BFIN, BOA, BVIN, C( 3) , 
:CANGL,CF(200) , CK, CLI , COF( 5 ) ,CW(200) , DEL (200) ,DMDOT(200) ,EA( 3,3) , 
: EPS (200) ,EZERO,F( 200,1) , FANGL , FNOBJ( 100 ) ,H(200) , HINF , HZ  (  200 )  , 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL( 100 ) ,R(200),RB(200), 
:RBASEI,R2I,RHOF(200) , ROOTI ( 4 ) ,RPM, ROOTR( 4 ) ,T(200) , TALFA, TB ( 20  0 ) , 
:TCC(200) ,TE(200) , THICK, THICKI , TIB ( 200 ) ,TINF , TS ( 200 ) ,TSAT,TSS, 
:TT(200) ,TZ,UF(200),X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,ZOA, 
:DOBF,DOTH,ICOR(200,3) , IFF, JINT, JLC, JTC, KFF( 50 ) ,KFIN( 50) , KT,NBAN, 
:NEL,NFIN,NSNP 

GUIDE  TO  FORTRAN  VARIABLE  NAMES 

AFOVAS  FIN  AREA/SMOOTH  AREA 

ALFA  FIN  HALF  ANGLE  (RADIANS) 

BFIN  HEIGHT  OF  FIN  (INCHES) 

BOA  TANGENT  OF  THE  FIN  HALF  ANGLE 

BVIN  HEIGHT  OF  FIN  (FEET) 

CALFA  COSINE  OF  ALFA 

CANGL  CONE  HALF  ANGLE  (DEGREES) 

CBASE  INSIDE  CIRCUMFERENCE  OF  CONDENSER  (FEET) 

CEXIT  INSIDE  CIRCUMFERENCE  AT  CONDENSER  EXIT  (FEET) 

CF  THERMAL  CONDUCTIVITY  OF  CONDENSATE  FILM  ( BTU/HR  FT  F) 

CL  CONDENSER  LENGTH  (FEET) 

CLI  CONDENSER  LENGTH  (INCHES) 

CPHI  COSINE  OF  PHI 

CRIT  CONVERGENCE  CRITERION 
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c 

DEL 

c 

DF 

c 

DIV 

c 

DMTOT 

c 

EPS 

c 

EPSEX 

c 

EPSO 

c 

EZERO 

c 

F 

c 

FANGL 

c 

G 

c 

H 

c 

HFG 

c 

IDG 

c 

IEL 

c 

IFF 

c 

IFIN 

c 

IFLUID 

c 

IGRAD 

c 

INFO 

c 

IONED 

c 

IOPT 

c 

IPRINT 

c 

ISTRAT 

c 

JINT 

c 

JLC 

c 

c 

JTC 

c 

c 

c 

c 

KFF 

c 

c 

KFIN 

c 

c 

c 

KT 

c 

NBAN 

c 

NBOTF 

c 

NBOTI 

c 

NCON 

c 

NDOBF 

c 

NDOTH 

c 

NDIV 

c 

NDV 

c 

NEFB 

c 

NEL 

c 

NEST 

c 

NRA 

c 

NRWK 

c 

NSNP 

c 

OBJ 

c 

OMEGA 

c 

PHI 

c 

PI 

c 

R2I 

c 

c 

RBASE 

c 

RBASEI 

c 

REXIT 

FILM  THICKNESS 

GRADIENT  OF  OBJECTIVE 

FLOATING  POINT  VALUE  OF  NDIV 

CONDENSATE  MASS  FLOW  RATE 

TROUGH  WIDTH  INCLUDING  INCREMENTAL  CHANGE 

TROUGH  WIDTH  AT  CONDENSER  EXIT 

TROUGH  WIDTH  AT  START  OF  CONDENSER 

FIN  BASE  WIDTH 

FORCE  VECTOR  OF  SYSTEM 

FIN  HALF  ANGLE  (DEGREES) 

CONSTRAINT  VALUES  ASSOCIATED  WITH  CURRENT  DESIGN 

CONVECTIVE  HEAT  TRANSFER  COEFFICIENT  ( BTU/HR  FT2  F) 

LATENT  HEAT  OF  VAPORIZATION  ( BTU/LBM ) 

CONSTRAINT  TYPE  IDENTIFIER 

THE  ELEMENT  NUMBER 

NO.  OF  ROWS  MINUS  ONE  OF  THE  UPPER  TRIANGULAR  FIN 

EQUALS  0  FOR  COPPER,  AND  1  FOR  STAINLESS  STEEL 

EQUALS  0  FOR  WATER,  AND  1  FOR  FREON 

GRADIENT  CALCULATION  IDENTIFIER 

CONTROL  PARAMETER 

ONE  DIMENSIONAL  SEARCH  IDENTIFIER 

OPTIMIZER  IDENTIFIER 

A  FOUR  DIGIT  PRINT  CONTROL 

OPTIMIZATION  STRATEGY  IDENTIFIER 

NO.  OF  COLUMNS  PLUS  ONE  BELOW  TRIANGULAR  FIN 

NUMBER  OF  SYSTEM  NODAL  POINT  LOCATED  AT 

THE  CENTER  OF  SYSTEM  COORDINATE 

NUMBER  OF  SYSTEM  NODAL  POINT  LOCATED  AT 

THE  JUNCTION  OF  THE  SYMMETRY  BOUNDARY  AND 

THE  LINE  OF  INTERSECTION  BETWEEN  THE  FIN 

AND  THE  CONDENSER  WALL 

NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  ALONG 

THE  FIN  CONVECTIVE  BOUNDARY 

NUMBER  OF  SYSTEM  NODAL  POINTS  LOCATED  ON  THE 

SYSMMETRIC  BOUNDARY  OF  TRIANGULAR  FIN  SECTION 

NOT  COUNTING  POINTS  AT  BASE  AND  APEX 

NUMBER  OF  COLUMNS  WITHIN  THE  TROUGH  WALL  SECTION 

SYSTEM  BAND  WIDTH 

LAST  ELEMENT  AT  BOTTOM  SIDE 

FIRST  ELEMENT  AT  BOTTOM  SIDE 

NUMBER  OF  CONSTRAINTS 

NUMBER  OF  ROWS  WITHIN  THE  FIN 

NUMBER  OF  ROWS  WITHIN  THE  TROUGH 

NUMBER  OF  INCREMENT 

NUMBER  OF  DESIGN  VARIABLES 

ELEMENT  NUMBER  AT  BASE  OF  FIN 

NUMBER  OF  ELEMENTS 

ELEMENT  NUMBER  AT  END  OF  TROUGH 

NUMBER  OF  ROWS  IN  ARRAY  A 

DIMENSIONAL  SIZE  OF  WK 

NUMBER  OF  SYSTEM  NODAL  POINTS 

VALUE  OF  THE  OBJECTIVE  FUNCTION  ASSOCIATED  WITH  X 

ROTATIONAL  SPEED  OF  HEAT  PIPE  ( RAD/HR ) 

CONE  HALF  ANGLE  (RADIANS) 

PI 

DISTANCE  FROM  CENTERLINE  OF  THE  HEAT  PIPE  TO  HALF  THE  F 

HEIGHT 

INSIDE  RADIUS  OF  CONDENSER  BASE  (FEET) 

INSIDE  RADIUS  OF  CONDENSER  BASE  (INCHES) 

INSIDE  RADIUS  OF  CONDENSER  EXIT  (FEET) 
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RPM  REVOLUTIONS  PER  MINUTE 

S  VECTOR  OF  DESIGN  VARIABLES 

SALFA  SINE  OF  ALFA 

SPHI  SINE  OF  PHI 

SURFAR  SURFACE  AREA 

THICK  CONDENSER  WALL  THICKNESS  (FEET) 

THICKI  CONDENSER  WALL  THICKNESS  (INCHES) 

TPHI  TANGENT  OF  PHI 

TZ  AMBIENT  TEMPERATURE 

UF  VISCOSITY 

VLB  LOWER  BOUNDS  ON  THE  DESIGN  VARIABLES  ' 

VUB  UPPER  BOUNDS  ON  THE  DESIGN  VARIABLES 

WK  REAL  WORK  ARRAY 

ZFIN  NUMBER  OF  FINS 

ZOA  RATIO  OF  TROUGH  WIDTH  TO  FIN  BASE  WIDTH 

ZSTAR  SURFACE  LENGTH  OF  THE  FIN  MINUS  THE  SURFACE  LENGTH 

COVERED  BY  THE  CONDENSATE  IN  THE  TROUGH 

ZZERO  SURFACE  LENGTH  OF  FIN 


PRINT*,  'INPUT  FILE  NAME' 
READ*,  NAME 
OPEN( 10 , FILE-NAME ) 
OPEN(15,FILE-'/HTPIPE  OUTPUT') 
OPEN( 14, FILE- '/DUMP  OUTPUT') 
OPEN( 13, FILE- '/GRAPH  OUTPUT') 


THE  FOLLOWING  READS  INPUT  DATA,  PERFORMS  HEAT 
TRANSFER  ANALYSIS,  AND  PRINTS  RESULTS. 


*****    INPUT  MODE    ***** 


ELEMENT  CONNECTIVITIES 

READ  (10,420)  NEL ,NSNP ,NBAN, IFLUID, IFIN 

WRITE  (15,430)  NEL ,NSNP , NBAN 

WRITE  (15,435)  IFLUID, IFIN 

WRITE  (15,436) 

WRITE  (15,437) 

READ  (10,440)  ( ICL , ( ICOR( IEL, I ) , 1-1 , 3 ) , IEL-1 ,NEL) 

WRITE  (15,450) 

THE  CONDENSER  GEOMETRY 

READ  (10,460)  CLI,CANGL,RBASEI,R2I, THICKI, BFIN,TZ 
WRITE  (15,470)  CLI,CANGL,RBASEI,R2I, THICKI, BFIN,TZ 
READ  (10,480)  NDIV, NEST,NEFB ,NBOTI , NBOTF 
WRITE  (15,490)  NDIV, NEST, NEFB , NBOTI , NBOTF 

DATA  FOR  RUNNING 

READ  (10,500)  RPM,TSS,TINF,HINF 
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WRITE  (15,510)  RPM,TSS,TINF,HINF 

C  THE  CONVERGENCE  CRITERIAN 

READ  (10,520)  CRIT 
WRITE  (15,530)  CRIT 

C  INTERNAL  FIN  GEOMETRY 

* 

READ  (10,54  0)  FANGL,NFIN 
WRITE  (15,550)  FANGL,NFIN 
WRITE(*,*)  FANGL,NFIN 
READ  (10,560)  IFF 
WRITE  (15,570)  IFF 

READ  (10,580)  ( KFIN( I ) , KFF( I ) , 1-1 , IFF ) 
READ  (10,590)  NDOBF,NDOTH, JTC, JLC, JINT, KT 
NHB-NEFB/2 
NBF-NBOTF+1 
DOBF  »  FLOAT(NDOBF) 
DOTH  »  FLOAT (NDOTH) 

WRITE  (15,600)  ICOR(NBOTI,2) , ICOR( NEFB, 1 ) , ICOR( NEST, 1 ) , 
1IC0R(NB0TF,1) 
* 

*  SET  CONSTRAINTS 
NRA-21 
NCOLA-30 
NRWK-  5000 
NRIWK-2000 
NDV-  1 
NCON-  2 
IGRAD-  0 

* 

*  INITIAL  DESIGN 
S(l)»  BFIN 

* 

*  BOUNDS 

VLB(l)-  0.0000001 
VUB(l)-  0.75 

* 

C      IDENTIFY  CONSTRAINTS 
C      NONLINEAR  CONSTRAINT 

IDG(1)-1 
C      LINEAR  CONSTRAINT 

IDG(2)=2 
* 

PRINT*,  'INPUT  THE  VALUES  FOR  ISTRAT, IOPT , IONED  AND  IPRINT' 

READ*,  ISTRAT, IOPT, IONED, IPRINT 
* 
* 
C      INITIALIZE  COUNTER 

NO=0.0 
C      CHANGE  THE  INTERNAL  PARAMETERS 

INFO— 2 

CALL  DADS ( INFO , ISTRAT , IOPT , IONED , I PRINT , IGRAD , NDV , NCON , S , VLB , 
:VUB,OBJ,G, IDG,NGT, IZ , DF ,W, NRA, NCOLA, WK , NRWK , IWK,NRIWK) 

IWK( 2)=0 

IWK(  3)=»200 

IWK(  5)=»4 

IWK(7)=500 

WK(3)=—  0.5 
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WK(6)«- 0.01 

WK(8)«0.05 

WK(9)»0.50 

WK(10)-0.05 

WK(ll)-0.005 

WK(13)«0.001 

WK(14)=0.0001 

WK(21)=0.004 

WK(22)=0.002 

WK(26)=0.004 

WK(37)=0. 0000000001 

)         CALL  DADS( INFO, ISTRAT , IOPT, IONED , IPRINT, IGRAD , NDV, NCON, S , VLB, 
:VUB,OBJ,G,IDG,NGT,IZ,DF,W,NRA,NCOLA,WK,NRWK,IWK,NRIWK) 
IF  (INFO.EQ.0)  GO  TO  360 

*****    EXECUTION  MODE    ***** 
NO-NO+1 

CONVERT  UNITS  OF  ALL  DIMENSIONAL  PARAMETERS 

TO  FEET.   CONVERT  UNITS  OF  ANGLES  TO  RADIANS. 

R2l-RBASEI-(S(l)/2) 

CL-CLI/12.0 

R2-R2I/12.0 

RBASE-RBASEI/12.0 

BVIN-S(1)/12.0 

DIV-FLOAT(NDIV) 

Pl'«3. 14159265358979 

PHI-2 . 0*CANGL*PI/360 . 0 

SPHI-SIN(PHI) 

CPHI-COS(PHI) 

TPHI-TAN(PHI) 

DELX-CL/DIV 

CBASE-2 . 0*PI *RBASE 

REXIT»RBASE+CL*TPHI 

CEXIT-2 . 0*PI*REXIT 

THICK-THICKI/12.0 

ALFA-FANGL*2 . 0*PI/360 . 0 

SALFA-SIN(ALFA) 

CALFA-COS(ALFA) 

TAL FA-TAN (ALFA) 

EZERO-2 . 0* ( S( 1 )/12 . 0 ) *TALFA 

ZOA=( (CBASE-(EZERO*NFIN) )/NFIN ) /EZERO 

BOUNDARY  CONDITIONS  AND  TEMPERATURE  ESTIMATES 
ALONG  THE  FIN  BOUNDARY 

DO  20  NTINF=NBOTI ,NBOTF 
20  TS(NTINF)=TINF 

DO  30  NNT-NBF,NEL 

TS(NNT)=0.0 
30  H(NNT)-0.0 

DO  40  IGT=1,NEST 

IE=*ICOR( IGTf 2) 
40  T(IE)=TZ 

IG=ICOR(NEST,l ) 

T( IG)=TZ 

OMEGA  IS  IN  RADIANS/HOUR 

OMEGA=RPM*2 . 0*PI*60  . 0 
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DO  50  KL=NBOTI,NBOTF 
50  H(KL)»HINF 

HIFN-HINF 

TSAT-TSS 

EPSO=ZOA*EZERO 

BOA-TALFA 

SURFAR-NFIN* (2.0*(S(l)/( 12*CALFA) )+EPSO) 

EPSEX»(CEXIT-(NFIN*EZERO) )/NFIN 

BETA«( EPSEX-EPSO)/DIV 

ZZERO=( S( 1 )/12 )/CALFA 

AFOVAS- ( ZOA+ ( 1 . /SALFA ) ) / ( 1 . +ZOA ) 

ZA-0,0 

DO  60  NSAT=1,NEST 
60  TS(NSAT)=TSAT 

TSOLID- ( TSAT+TINF )/2 . 0 
C        TEMPORARY  CHANGE  -  TFILM 

ASMOOTH-0.0 

ACASE-0.0 

QT-0.0 

OBJ-0.0 

QT1-0.0 

QTF-0.0 

QTRF-0.0 

QTOT-0.0 

DMTOT-0.0 

NK-NDIV+1 

DO  350  NI»1,NK 
C      R  IS  THE  INCREMENTAL  CHANGE  IN  THE  RADIUS  OF  THE  CONDENSER 

R(NI)«R2+NI*DELX*SPHI 

RB ( NI ) -RBASE+NI *DELX*SPHI 
C      EPS  IS  THE  INCREMENTAL  CHANGE  IN  THE  TROUGH  WIDTH 

EPS ( NI ) -EPSO+NI *BETA 

APS-EPS(NI) 
C 

C  NODAL  POINT  COORDINATES 

C 

CALL  COORD 
65   Z(1)-ZA 

DO  70  IZEL«1,NEFB 

NA»ICOR(IZEL,l) 

NB*ICOR(IZEL,2) 

XE«X(NA)-X(NB) 

YE=»Y(NA)-Y(NB) 

ELZ=»SQRT(XE**2+YE**2) 
70  Z(IZEL+1)=Z(  IZED+ELZ 

XZB=X(ICOR(NHB,l) ) -X( ICOR( 1 , 2 ) ) 

YZB=»Y(ICOR(NHB,l) )-Y( ICOR(l,2) ) 

ZB=SQRT(XZB**2+YZB**2) 

ZC-ZZERO 

IM-1 
C 

C        PARABOLIC  TEMPERATURE  DISTRIBUTION  ALONG  THE  FIN 
C        BOUNDARY, USING  LAGRANGE  INTERPOLATION 
C 

80  TPl=>T( ICOR(l, 2) ) 

TP2=T( ICOR(NHB,l ) ) 

TP3=T( ICOR(NEFB,l) ) 

APl=TPl/(ZB*ZC) 

AP2=TP2/(ZB*(ZB-ZC) ) 

AP3=TP3/(ZC*(ZC-ZB) ) 
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BPl— (ZB+ZC)*APl 

BP2  —  ZC*AP2 

BP3  — ZB*AP3 

A1-AP1+AP2+AP3 

B1-BP1+BP2+BP3 

TC-0.0 

DO  90  NY=1,NEST 

90  TOTC+T(ICOR(NY,2)  ) 
AY-FLOAT(NY+l) 

TF-(TC+T( I COR ( NY, 1) )+AY*TS(NY) )/(2.0*AY) 

SOLID-FLUID  PROPERTIES 
WATER  PROPERTIES 

IF(IFLUID.EQ.l)  GO  TO  91 

HFG-1093.88-0.5703*TS(1)+0.00012819*(TS(1)**2) 
1-0. 0000008824 *(TS(1)**3) 

RHOF(NI ) -62. 774-0. 00255698 *TF-0. 000053 572*TF**2 

CF(NI)-0. 3034+0. 0 007 38927 *TF-0. 000001 47321 *TF**2 

UF(NI)-0.001397-0.000014669*TF+0.00000  006  31253*TF**2-0.000  00 
100000976569*TF**3 

FREON  PROPERTIES 

91  IF(IFLUID.EQ.O)  GO  TO  92 

HFG-69. 5459-0. 0156011*TS( 1 )-0 . 000455294* (TS( 1 )**2)+0 . 00000104144* ( 
1TS(1)**3) 

RHOF (NI) -102. 059-0. 025364 *TF-0. 000 502649* (TF**2) +0.00 0001 3 54 07 * ( TF 
1**3) 

CF(NI )-0. 0  594858-0. 00042976 5*TF+0 . 00000 3 48218*TF**2-0 . 000000010416 
18*TF**3 

UF(NI) -0.0 0078-0. 00000 52 5 *TF+0. 0000000125 *TF**2 

92  UF(NI)-3600*UF(NI) 

IF( IFIN.EQ.l)  GO  TO  93 

CW(NI) -2 3 1.7772-0. 02222 *TSOLID 

93  IF(IFIN.EQ.O)  GO  TO  94 
CW(NI)-8.776+0.00265*TSOLID 

94  CK-CW(NI) 

CONST-RHOF(NI ) **2*OMEGA**2*HFG*CPHI *CALFA*R( NI ) 

INITIAL  FILM  THICKNESS 
IF  (NI.GT.l)  GO  TO  100 

DEL(1)-1.107*( ( (TSAT-TINF)*CF(NI)/(UF(NI)*HFG) ) ** . 25 ) * ( ( UF ( NI ) /( 
1RH0F(NI )*OMEGA) )**0.5) 
100  CONTINUE 

AVERAGE  ELEMENT  CONVECTIVE  COEFFICIENT  ALONG 
THE  FIN  BOUNDARY 

ZSTAR=ZZERO-DEL(NI )/CALFA 

AZZ-DEL(NI)/SALFA 

ZZ-ZSTAR 

HDEN=( (-Al*ZZ**3)/3.0-(Bl*ZZ**2)/2.0)+ZZ*(TSAT-T(l) ) 

AZS=ABS( 4*CF(NI )*UF(NI ) *HDEN/CONST ) * *0 . 2 5 

HAC=0.0 

DO  190  IEL=1,NEFB 

AZ=Z( IEL) 

BZ=Z( IEL+1) 

IF    ( ZSTAR.LE.BZ )    GO    TO    110 
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GO  TO  120 
110  IF  (HAC.NE.0.0)  GO  TO  18  0 

BZ-ZSTAR 
120  IF  (IEL.NE.l)  GO  TO  130 

AK=-(BZ-AZ)/5.0 

ZZ-AK 

GO  TO  14  0 
130  AK-(BZ-AZ)/4.0 

ZZ-AZ 
140  ZEL=4*AK 

DO  150  NH=1,5 

HDEN=(-1.0*(Al*ZZ**3/3.0+Bl*ZZ**2/2.0) ) +ZZ* ( TSAT-T( 1 ) 

HZ(NH)»ABS(CF(NI)**3*CONST/(4*UF(NI)*HDEN) )**0.25 
150  ZZ-ZZ+AK 
C      AVERAGE  H  USING  SIMPSONS  RULE 

CONH»AK*(HZ(l)+4*HZ(2)+2*HZ(3)+4*HZ(4)+HZ(5) )/(3*ZEL) 

IF  (ZSTAR.EQ.BZ)  GO  TO  160 

H(IEL)-CONH 

GO  TO  190 
160  AZ-ZSTAR 

HAZ-CONH*(AZ-Z(IEL) ) 

DELA-AZS 
170  BZ-Z(IEL+1) 

DELB-(BZ-ZSTAR)*AZZ/(ZZERO-ZSTAR) 

DELZ- ( DELA+DELB ) /2 . 0 

HAC« ( BZ-AZ ) *CF ( NI ) /DELZ 

H(IEL)-(HAZ+HAC)/(BZ-Z(IEL) ) 

GO  TO  190 
180  AZ-Z(IEL) 

DELA-DELB 

HAZ-0.0 

GO  TO  170 
190  CONTINUE 

NETI-NEFB+1 

DO  200  IEL-NETI,NEST 
200  H(IEL)=CF(NI)/DEL(NI) 
C 

C  ENTRY  INTO  THE  FINITE  ELEMENT  SOLUTION 

C 

CALL  FORMAF 

CALL  BANDEC  ( NSNP , NBAN , 1 ) 
C 

C  THE  TEMPERATURE  DISTRIBUTION 

C 

DO  210  NT=1,NSNP 
210  T(NT)»F(NT,1) 

TIB(NI)=»T(ICOR(NBOTIf  2)  ) 

TT(NI)=T( ICOR(NEFB,l) ) 

TE(NI)=T(ICOR(NEST,l) ) 

TB(NI)«T(ICOR(NBOTF,l) ) 

TTS-0.0 

DO  220  NS=»1,NSNP 
220  TTS-TTS+T(NS) 

PN-FLOAT(NS) 

TSOLID=TTS/PN 
C 

C  Q  AT  THE  BOTTOM  SIDE 

C 

QBI=0.0 

DO  230  IBEL=NBOTI/NBOTF 
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NKA-IC0R(IBEL,1) 
NKB-ICOR( IBEL,2) 
XB-X(NKA)-X(NKB) 
YB-Y(NKA)-Y(NKB) 
ELB=SQRT(XB**2+YB**2) 
230  QBI*QBI+(T(NKA)+T(NKB)-2*TS( I BEL) )*ELB*H( IBEL)/2.0 
QB(NI  )=-QBI*DELX 

UNTIL  CONVERGENCE  CRITERIA  IS  MET 

TO  2  40 


ITERATION 

IF 

(IM 

.  EQ. 

1) 

GO 

QJ* 

■QBI 

GO 

to  : 

250 

QI» 

■QBI 

240 

IM-2 

GO  TO  80 
250  AQ-ABS(QJ-QI)/QJ 

IF  (AQ.LE.CRIT)  GO  TO  260 

QI-QJ 

GO  TO  80 
260  DMDOT(NI)=»2.*QBI*DELX/HFG 

DMTOT«DMTOT+DMDOT( NI ) 

Cl-RHOF(NI)**2*OMEGA**2*R(NI)*SPHI/(3*UF(NI) ) 

XCOF(l)»-DMTOT 

XCOF(2)-0.0 

XCOF(3)-0.0 

XCOF(4)-Cl*EPS(NI) 

XCOF(5)-Cl*TALFA 

M-4 

CALL  DPOLRT  (M,IER) 

IF  (ROOTR(l) .GT.0.0)  GO  TO  270 

IF  (ROOTR(2) .GT.0.0)  GO  TO  280 

IF  (ROOTR( 3) .GT.0.0)  GO  TO  290 

IF  (ROOTR( 4) .GT.0.0)  GO  TO  300 

THE  CONDENSATE  THICKNESS 

270  DEL(NI+l)-ROOTR( 1) 

GO  TO  310 
280  DEL(NI+l)»ROOTR(2) 

GO  TO  310 
290  DEL(NI+l)»ROOTR( 3) 

GO  TO  310 
•300  DEL(NI  +  l)=»ROOTR(4) 
310  QEL-0.0 

IF  (NI.NE.l)  GO  TO  320 

Q  FROM  THE  TOP  SIDE 
Q  THROUGH  FIN 

QEL-0.0 
320  DO  330  IQEL«1,NEFB 

KA»ICOR(IQEL,l) 

KB=ICOR(IQEL, 2) 

XQEL-X(KA)-X(KB) 

YQEL=Y(KB)-Y(KA) 

ELM»SQRT(XQEL**2+YQEL**2) 

QEL-QEL+(2*TS(  IQEL  ) -T(  KA) -T(  KB  )  )*ELM*H(  IQED/2.0 
330    CONTINUE 
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QINC(NI)-QEL*DELX 

AMTOT(NI)»DMTOT 

QET-QEL*DELX*NFIN*2 

QT-QT+QET 

QA»QBI*DELX*NFIN*2 
C 

C      Q  THROUGH  TROUGH 
C 

QTRF-0.0 

DO  340  IQEL=NEFB+1,NEST 

KA-ICOR(IQEL,l) 

KB»ICOR(IQEL,2) 

XQEL«X(KA)-X(KB) 

YQEL-Y(KB)-Y(KA) 

ELM-SQRT(XQEL**2+YQEL**2) 

QTRF-QTRF+ ( 2  *TS ( IQEL ) -T ( KA ) -T ( KB ) ) *ELM*H ( IQEL ) /2 . 0 
340  CONTINUE 

QTINC ( NI ) -QTRF*DELX 

QTOTAL ( NI ) -QINC ( NI ) +QTINC ( NI ) 

QTRFT=QTRF*DELX*NFIN*2 . 

QTF-QTF+QTRFT 

QTOT-QTOT+QA 

ASMOOTH«2*PI*RB(NI)*DELX+ASMOOTH 

ACASE-ACASE+NFIN*( ( (2*S(1) )/(12*CALFA) )+EPSO)*DELX 
350  CONTINUE 
C      EVALUATE  OBJECTIVE  FUNCTIONS  AND  CONSTRAINTS 

OBJ-(QTOT) 

WRITE(15,*)  'OBJECTIVE  FUNCTION* ',    OBJ^BFIN-',  S(l) 
C      FIRST  CONSTRAINT  IS  TO  ENSURE  THE  RATIO  ZOA  IS  NOT  NEGATIVE 

6(1)— ( ( (CBASE-(2*NFIN*S(1)*TALFA/12) )/NFIN)/( S ( 1 ) *TALFA/6 ) ) 

G(1)-1000.0*G(1) 
C      THE  SECOND  CONSTRAINT  IS  TO  ENSURE  THE  CONDENSATE  LEVEL  IS  NO 
C      GREATER  THAN  THE  FIN  HEIGHT 

G(2)  — (  (S(1)/12.0)-DEL(NI)) 

XPLOT(NO)-Sd) 

FNOBJ(NO)— OBJ 

ARATIO-ACASE/ASMOOTH 

WRITE (13,*)  XPLOT(NO) ,FNOBJ(NO) 

WRITE (14,*)  ARATIO,  FNOBJ(NO) 

GO  TO  10 
360  CONTINUE 
C 

BFIN-SU) 
C      *****    OUTPUT  MODE    ***** 

WRITE  (15,630) 

DO  370  NR=»1,NK 
370  WRITE  (15,640)  NR,QINC( NR ), QTINC ( NR) , QTOTAL ( NR) 

WRITE  (15,650)  QT,QTF 

WRITE  (15,660)  CLI , CANGL, RBASEI , R2I , THICKI , BFIN, RPM, TSS , TINF , 
1RIT,FANGL,Z0A,IFF 

WRITE  (15,661)  AFOVAS 

WRITE  (15,670)  BOA,ZOA,NFIN,BVIN,SURFAR 

WRITE  (15,680) 

DO  380  NP»1,NSNP 

TCC (  NP  )=..  5555555 *(T(NP)  -32) 
380  WRITE  (15,690)  NP , X( NP ) , Y( NP ) , T( NP ) , TCC( NP ) 

WRITE  (15,700) 

DO    390    KKL-l,NBOTF 

NKX-ICOR(KKL,l) 

NKY=ICOR(KKL,2) 
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XP-X(NKX)-X(NKY) 

YP-Y(NKX)-Y(NKY) 

EXY-SQRT(XP**2+YP**2) 

QEP-ABS( (T(NKX)+T(NKY)-2*TS(KKL) ) *EXY*H( KKL)/2 . 0 ) 

QEP«QEP*DELX 

KKL,H(KKL) ,EXY,QEP 

CRIT 

HFG , NFIN  f  H ( NBOTF )  , TSAT , RPM , QTOT , QT , FANGL 


390  WRITE  (15,710) 

WRITE  (15,720) 

WRITE  (15,730) 

WRITE(15,734) 

WRITE(15,735) 

WRITE(15,735) 

WRITE( 15,735) 

WRITE(15,735) 

WRITE  (15,740) 

DO  400  NR-1,NK 
400  WRITE  (15,750) 
1TB (NR) 

WRITE  (15,760) 

DO  410  NG=1,NDIV,2 
410  WRTTE  (15,770)  NG, CW( NG) , CF( NG 
lQINC(NG) 

RETURN 


ROOTR(l) ,ROOTI(l) 
ROOTR(2) ,ROOTI(2) 
ROOTR(3) ,ROOTI(3) 
ROOTR(4) ,ROOTI(4) 


NR,DEL(NR) , QB ( NR ) ,AMTOT(NR) ,TIB(NR) ,TT(NR) ,TE(NR) 


, RHOF ( NG ) , UF ( NG ) , EPS ( NG ) , R ( NG 


412 
420 
430 

435 
436 

437 
440 
450 
460 
470 


FORMAT 
FORMAT 
FORMAT 
115,/, 2X 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


15,/,4X,5HR2I-  ,E12.5,/,4X,8HTHICKI-  , E12 . 5 ,/, 4X, 6HBFIN=  ,El2.5,/,4 


480 
490 

500 
510 

520 
530 
540 
550 
560 
570 
580 
590 
600 

610 
620 
630 
640 
650 
660 


2X,4HTZ» 

FORMAT 

FORMAT 

1HNBOTI' 
FORMAT 
FORMAT 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


1,6X,3HTB-,I5) 


8X,E12.5,8X,E12.5) 

515) 

/2X , 1 5HNO . OF . ELEMENTS- , I 5 , 1 OX , 1 8HNO . OF . SYSTEM  N . P . = , 

13HNO.OF  BANDED-, 15) 

/2X, 'I FLUID-' ,I5,10X, ' IFIN-' ,15) 

2X,' I FLUID  -  0  FOR  WATER,  AND  1  FOR  FREON' ) 

2X,'IFIN  -  0  FOR  COPPER,  AND  1  FOR  STAINLESS  STEEL') 

415) 

/2X, 7HELEMENT, 10X, 3HNP1 , 14X, 3HNP2 , 15X, 3HNP3 ) 

7G10.5) 

4X,5HCLI-  ,E12.5,/,4X,7HCANGL=«  , E12 . 5 ,/, 4X , 8HRBASEI . = , E12 


,E12.5) 

515) 

4X,6HNDIV-  ,I10,/,4X,6HNEST-  , 110 ,/, 4X, 6HNEFB=»  ,I10,/,4X,7 

,I10,/,4X,7HNBOTF-  ,110) 

4F10.2) 

4X,5HRPM»  ,E12.5,/,4X,5HTSS=  , E12 . 5 ,/, 4X, 6HTINF=  ,E12.5,/, 


14X,6HHINF-  ,E12.5) 


G10.9) 

4X,6HCRIT=  ,E12.5) 

2G10.5) 

4X,7HFANGL=  , E12 . 5 ,/, 4X , 6HNFIN=  ,15) 

15) 

4X,5HIFF=  ,110) 

1615) 

615) 

///5X,4HTIB«,I5,10X,3HTT-,I5,/,5X,3HTE-,I5,/ 


FORMAT  ( //l OX, 17HCRASH, CRASH, CRASH) 

FORMAT  (//5X,4(E12.7,3X) ) 

FORMAT  ( 2X , 7HSTATION , 2X , 4HQFIN , 17X , 7HQTROUGH , 1 5X , 6HQTOTAL ) 

FORMAT  (4X, 1 5, El 2. 5,1  OX, El 2. 5,1  OX, El 2. 5) 

FORMAT  (//,4X,11HQFIN  TOTAL= , E12 . 5 , 10X , 1 5HQTROUGH  TOTAL=  , E12.5) 

FORMAT  (/////, 4X,5HCL I =  , E12 . 5 , 5X , 7HCANGL-  , E12 . 5 , /, 4X, 8HRBASEI=  , 
1E12.5,2X,5HR2I=  , E12 . 5 ,/, 4X, 8HTHICKI=  , E12 . 5 , 2X , 6HBFIN=  ,El2.5,/,4 
2X,5HRPM-  ,E12.5,5X,5HTSS-  , E12 . 5 ,/, 4X, 6HTINF=  , E12 . 5 , 4X, 6HHINF=  ,E 
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312.5,/, 4X,6HCRIT=  , El2 . 5 , 4X, 7HFANGL-  , E12 . 5 ,/, 4X, 5HZ0A-  ,E12 

45HIFF-  ,110) 
661  FORMAT (4X, 'FIN  AREA/SMOOTH  AREA- ' , El 2 . 5 ) 
670  FORMAT  (1H1,//2X,4HB0A-,G12.5,5X,4HZ0A=,G12.5,5X,5HNFIN«,I5, 

1/, 5HBVIN- , G12 . 5 , 5X, 13HSURFACE  AREA- , G12 . 5 ) 
680  FORMAT  (//5X, 2HNP , 6X, 1HX, 12X, 1HY, 12X, lHT, 12X, 2HTC ) 
690  FORMAT  (/2X,I3,3X,4(F10.6,3X) ) 

700  FORMAT  ( //2X, 2HEL , 8X, 1HH , 11X, 9HEL-LENGTH , 15X, 4HQ-EL) 
710  FORMAT  (/2X, 12 f 3Xf E12 . 5 , 3X, E12 . 5 , 10X, E12 . 5 ) 
720  FORMAT  ( /2X, 22HCONVERGENCE  CRITERIAN- , E15 . 8 ) 
730  FORMAT  ( 1H  ,//, 5X, 4HHFG= , E12 . 5 ,/, 5X, 11HNO . OF  FINS-, 15 ,/, 5X, 

1 6HH-OUT- , El 2 . 5 , / , 5X , 5HTSAT- , El 2 . 5 , / , 5X , 4HRPM- , El 2 . 5 , / , 5X , 6HQ 

2E12.5,/,5X,6HQFIN  - , El 2 . 5 ,/, 5X, 11HHALF-ANGLE- , F8 . 3 ) 

734  FORMAT (/5X, 'ROOTS: ' ,5X, 'REAL  PARTS ', 15X, ' IMAGINARY  PARTS') 

735  F0RMAT(15X,E12.5,15X,E12.5) 

740  FORMAT  ( 1H0 , 6X, 1HJ, 4X, 14HFILM  THICKNESS , 8X, 2HQB , 10X, 8HMASS-T 

1/, 4X, 3HTIB, 8X, 2HTT, 10X, 2HTE , 8X, 2HTB ) 
750  FORMAT  (1H  , 4X, 14 , 4X, F12 . 10 , 4X, F10 . 4 , 6X, F9 . 5 , 6X, F5 . 1 ,/, 

16X,F5.1,6X,F5.1,6X,F5.1) 
760  FORMAT  (1H0,6X,1HJ,6X,6HK-WALL,4X,6HK-FILM,3X,7HDENSITY,4X,9 

1FILM,/,6X,7HEPSIL0N,5X,6HRADIUS,5X,4HQINC) 
770  FORMAT  ( 1H  , 4X, 14 , 4X, F7 . 3 , 4X, F6 . 4 , 4X, F6 . 3 , 4X, F9 . 7 , 

1/, 4X , F9 . 7 , 4X , F7 . 5 , 5X , F5 . 1 , IX , F7 . 3 ) 
END 


* 


*      THIS  SUBROUTINE  DEFINES  THE  POSITIONS  OF  SYSTEM  COORDINATE  P 

* 

SUBROUTINE  COORD 

* 

COMMON/ADS/DF(21) ,G(10) ,IDG(100) , IGRAD, INFO, IOPT, IONED, IPRIN 
:ISTRAT,IWK(2000) ,IZ(30) , OB J , S ( 2 ) ,VLB(2) ,VUB(2) ,W( 21 , 30 ) ,WK( 5 
: NCOLA , NCON , NDV , NGT , NRA , NR IWK , NRWK 

* 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B( 3) , BFIN, BOA, BVIN, C( 3 
:CANGL,CF(200) , CK, CLI , COF( 5 ) ,CW(200) , DEL (200) ,DMDOT(200) ,EA(3 
: EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ( 100 ) , H( 2  00 ) , HINF , HZ ( 200 ) 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL( 100 ) ,R(200),RB(200) 
:RBASEI,R2I,RHOF(200) ,ROOTI(4) , RPM, ROOTR( 4 ) ,T(200) ,TALFA,TB(2 
:TCC(200) ,TE(200) , THICK, THICKI , TIB ( 200 ) , TINF , TS ( 200 ) ,TSAT,TSS 
:TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200),Z 
:DOBF,DOTH,ICOR(200,3) , IFF, JINT, JLC , JTC, KFF( 50) ,KFIN( 50) ,KT,N 
:NEL,NFIN,NSNP 
* 
* 
C      DELH  IS  THE  STANDARD  DIVISION  OF  FIN  HEIGHT 

DELH-S(l)/(12*DOBF) 

X(1)=0.0 

Y(1)-THICK+(S(1)/12) 

N-l 

DO  20  I-1,IFF 

ICA-KFIN( I) 

ICB-KFF( I) 

CBA- FLOAT ( ICB-ICA) 

AN-0.0 

DO  10  II-ICA,ICB 

X( II )=X( 1 )+N*AN*DELH*TALFA/CBA 

Y(II)=Y(1)-N*DELH 


10  AN-AN+1.0 
20  N=N+1 
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AN-0.0 

ICD-ICB-ICA+1 

DO  50  J-JTC, JLC, JINT 

X(J)-X(1) 

Y( J)-(1.0-AN/DOTH)*THICK 

DO  30  JJ-1,ICD 

X(  J+JJ ) «X( J ) +JJ*EZERO/( 2* ( CBA+1 . 0 )  ) 
30  Y(  J+JJ)=»Y(  J) 

JJ-ICD 

DO  40  K-1,KT 

X( J+JJ+K)«X( J+JJ) +K*APS/( 2.0* KT) 
40  Y( J+JJ+K)-Y( J) 
50  AN-AN+1.0 

RETURN 

END 


THIS  SUBROUTINUE  IS  USED  TO  FORMULATE  THE  FEM  EQUATIONS 

SUBROUTINE  FORMAF 

COMMON/ADS/DF(21) ,G(10) , IDG(100) , IGRAD, INFO, IOPT, IONED , IPRINT, 
:ISTRAT,IWK(2000) , IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W(21,30),WK(5000) 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B(3) , BFIN, BOA, BVIN, C( 3) , 
:CANGL,CF(200) , CK, CLI , COF( 5 ) ,CW(200) , DEL (200) ,DMDOT(200) ,EA(3,3) , 
: EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ( 100 ) ,H(200) ,HINF, HZ ( 200 )  ,' 
:QB(200) ,QINC(200) ,QTINC(200) , QTOT , QTOTAL ( 100 ) ,R( 200 ) , RB ( 200 ) , 
:RBASEI,R2I,RHOF(200) ,ROOTI(4) , RPM,ROOTR( 4 ) ,T(200) , TALFA, TB ( 200 ) , 
:TCC(200) ,TE(200) , THICK, THICKI , TIB ( 200 ) , TINF, TS ( 200 ) ,TSAT,TSS,. 
:TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,  ZOA, 
:DOBF,DOTH,ICOR(200,3) , IFF, JINT, JLC, JTC, KFF( 50 ) ,KFIN( 50) , KT,NBAN, 
:NEL,NFIN,NSNP 

DO  20  N-1,NSNP 

F(N,l)-0.0 

DO  10  MA»1,NBAN 
10  A(N,MA)-0.0 
20  CONTINUE 

DO  70  IEL-1,NEL 

IA»ICOR(IEL,l) 

IB-ICOR(IEL,2) 

IC-ICOR( IEL,3) 

B(1)-Y(IB)-Y( IC) 

B(2)-Y( IC)-Y(IA) 

B( 3)-Y( IA)-Y(IB) 

C(1)-X(IC)-X(IB) 

C(2)-X(IA)-X(IC) 

C(3)«X(IB)-X(IA) 

LENGTH  BETWEEN  ELEMENT  NODES  1  AND  2 

EL-SQRT(C( 3)**2+B( 3)**2) 

AREA  OF  TRIANGULAR  ELEMENT 

AS=ABS( (B(1)*C(2)-B(2)*C(1) )/2.0) 

HC=H( IEL)/CK 

DO    60    J=»l,3 

JJ=»ICOR( IEL, J) 

DO    50    K=l,3 

KK-ICOR(IEL,K) 
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FORMING  THE  A  MATRIX 

EA(J,K)-(B(J)*B(K)+C(J)*C(K) )/(4*AS) 

IF  (HC.EQ.0.0)  GO  TO  40 

HEL=HC*EL/6.0 

IF  (J.EQ.3)  GO  TO  4  0 

IF  (K.EQ.3)  GO  TO  4  0 

IF  (J.EQ.K)  GO  TO  3  0 

EA( J,K)=EA( J,K)+HEL 

GO  TO  4  0 
30  EA( J,K)=EA( J,K)+2*HEL 
40  IF  (KK.LT.JJ)  GO  TO  50 

NW-KK-JJ+1 

A( JJ,NW)»A( JJ,NW)+EA( J,K) 
50  CONTINUE 
60  CONTINUE 

FORMING  THE  F  MATRIX 

FE-HC*TS ( IEL) *EL/2 . 0 

F(IA,1)=F(IA,1)+FE 

F(IB,1)=F(IB,1)+FE 
70  CONTINUE 

RETURN 

END 


EQUATION  SOLVER  OF  A  SYMMETRIC  MATRIX  THAT  HAS  BEEN  TRANS- 
FORMED INTO  BANDED  FORM. 

SUBROUTINE  BANDEC  ( NEQ, MAXB,NVEC) 

COMMON/ADS/DF(21) ,G(10) , IDG(100) , IGRAD, INFO, IOPT, IONED, IPRIN 
ISTRAT,IWK(2000) ,IZ(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W( 21 , 30 ) ,WK ( 5 
NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B( 3) , BFIN, BOA, BVIN, C( 3 
CANGL,CF(200) ,CK,CLI,COF( 5) ,CW(200) , DEL (200) ,DMDOT(200) ,EA( 3 
EPS (200) ,EZERO,F( 200,1) , FANGL , FNOB J ( 100 ) ,H(200) ,HINF,HZ( 200) 
QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL( 100 ) ,R(200),RB(200) 
RBASEI,R2I,RHOF(200) ,ROOTI(4) , RPM, ROOTR( 4 ) ,T(200) ,TALFA,TB(2 
TCC(200) ,TE(200) , THICK, THICKI , TIB( 200 ) , TINF, TS ( 200 ) ,TSAT,TSS 
TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200),Z 
DOBF,DOTH,ICOR(200,3) , IFF , JINT, JLC , JTC , KFF( 50 ) , KFIN( 50 ) , KT , N 
NEL,NFIN,NSNP 


LOOP=NEQ-l 

DO  20  1=1, LOOP 

MB-I+1 

NB-MIN0 ( I+MAXB-1 , NEQ ) 

DO  20  J=MB,NB 

L-J+2-MB 

D=A(I,L)/A(I,1) 

DO  10  MM-1,NVEC 
10  F( J,MM)-F( J,MM)-D*F( I , MM ) 

MM-MIN0 ( MAXB-L+1 , NEQ- J+l ) 

DO  20  K=1,MM 

NN-L+K-1 
20  A( J,K)=A( J,K)-D*A( I , NN ) 

DO  30  I=1,NVEC 
30  F(NEQ,I)=F(NEQ,I)/A(NEQ,1) 

DO  50  1=2, NEQ 
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J-NEQ-I+1 

K-MINO ( NEQ- J+l , MAXB ) 

DO  50  MM-1,NVEC 

DO  40  L»2,K 

MB-J+L-1 
4  0  F ( J , MM ) =  F ( J , MM ) -A ( J , L) *  F ( MB , MM ) 
50  F( J,MM)-F( J,MM)/A( J , 1 ) 

RETURN 

END 

SUBROUTINE  DPOLRT  COMPUTES  THE  ROOTS  OF  A  REAL 
POLYNOMIAL  USING  A  NEWTON-RAPHSON  ITERATIVE 
TECHNIQUE. 

SUBROUTINE  DPOLRT  (M,IER) 

COMMON/ADS/DF(21) ,G(10) , IDG(IOO) , IGRAD , INFO , IOPT , IONED , IPRINT , 
:ISTRAT,IWK( 2000) ,12(30) ,OBJ,S(2) ,VLB(2) ,VUB(2) ,W(21,30),WK(5000) 
: NCOLA , NCON , NDV , NGT , NRA , NRIWK , NRWK 

COMMON/OLLIE/A(200,50) ,AMTOT(200) ,APS,B(3) , BFIN, BOA, BVIN, C( 3 ) , 
:CANGL,CF(200) , CK, CLI , COF( 5 ) ,CW(200) , DEL (200) ,DMDOT(200) ,EA( 3,3) , 
: EPS (200) ,EZERO,F( 200,1) , FANGL, FNOBJ( 100 ) ,H(200) , HINF,HZ ( 200 ) , 
:QB(200) ,QINC(200) ,QTINC(200) ,QTOT,QTOTAL( 100 ) ,R(200),RB(200), 
:RBASEI,R2I,RHOF(200) ,ROOTI(4) ,RPM,ROOTR( 4 ) ,T(200) ,TALFA, TB( 200 ) , 
:TCC(200)  ,TE(200)  , THICK, THICKI ,  TIB (  200.)  ,  TINF,TS(  200  )  ,TSAT,TSS, 
:TT(200) ,TZ,UF(200) ,X(200) ,XCOF(5) ,XPLOT(200) ,Y(200),Z(200) ,ZOA, 
: DOBF , DOTH , ICOR ( 20  0 , 3 ) , I FF , JINT , JLC , JTC , KFF ( 5  0 ) , KFIN ( 5  0 ) , KT , NBAN , 
:NEL,NFIN,NSNP 

IFIT-0 
N-M 
IER-0 

IF  (XCOF(N+l))  10,40,10 
10  IF  (N)  20,20,60 

SET  ERROR  CODE  TO  1 

20  IER-1 
30  RETURN 

SET  ERROR  CODE  TO  4 

40  IER-4 

GO  TO  30 

SET  ERROR  CODE  TO  2 

50  IER-2 

GO  TO  3  0 
60  IF  (N-36)  70,70,50 
70  NX=N 

NXX-N+1 

N2-1 

KJl-N+1 

DO  80  L-1,KJ1 

MT=-KJ1-L+1 
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80  COF(MT)-XCOF(L) 
C 

C  SET  INITIAL  VALUES 

C 

90  XO-. 00500101 
YO-0. 01000101 
C 

C  ZERO  INITIAL  VALUE  COUNTER 

C 

IN-0 
100  XX=XO 
C 

C         INCREMENT  INITIAL  VALUES  AND  COUNTER 
C 


XO—  10.0*YO 

YO— 10.0*XX 
C 

C         SET  X  AND  Y  TO  CURRENT  VALUE 
C 

XX=  XO 

YY-YO 

IN-IN+1 

GO  TO  120 
110  IFIT-1 

XPR»XX 

YPR-YY 
C 

C         EVALUATE  POLYNOMIAL  AND  DERIVATIVES 
C 

120  ICT-0 
130  UX-0.0  . 

UY-0o0 

V-0.0 

YT»0.0 

XT-1.0 

U-COF(N+l) 

IF  (U)  140,270,140 
140  DO  150  I-1,N 

L-N-I+l 

XT2-XX*XT-YY*YT 

YT2»XX*YT+YY*XT 

U-U+COF(L)*XT2 

V=V+COF(L)*YT2 

FI  =  I 

UX=UX+FI*XT*COF( L) 

UY=UY-FI*YT*COF(L) 

XT-XT2 
150  YT»YT2 

SUMSQ=UX*UX+UY*UY 

IF  (SUMSQ)  160,230,160 
160  DX»(V*UY-U*UX)/SUMSQ 

XX-XX+DX 

DY—  (  U*UY+V*UX  )  /SUMSQ 

YY-YY+DY 

IF  (ABS(DY)+ABS(DX)-1.0E-05)  210,170,170 
C 

C  STEP  ITERATION  COUNTER 

C 

170  ICT-ICT+1 

IF  (ICT-500)  130,180,180 
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180  IF  (IFIT)  210,190,210 
190  IF  (IN-5)  100,200,200 

SET  ERROR  CODE  TO  3 

200  IER-3 

GO  TO  30 
210  DO  220  L-1,NXX 

MT-KJl-L+1 

TEMP-XCOF(MT) 

XCOF(MT)-COF(L) 
220  COF(L)»TEMP 

ITEMP-N 

N-NX 

NX-ITEMP 

IF  (IFIT)  250,110,250 
230  IF  (IFIT)  240,100,240 
240  XX-XPR 

YY-YPR 
250  IFIT-0 

IF  (ABS(YY/XX)-1.0E-04)  280,260,260 
260  ALPHA- XX+XX 

SUMSQ-XX*XX+YY*YY 

N-N-2 

GO  TO  290 
270  XX-0.0 

NX-NX- 1 

NXX-NXX-1 
280  YY-0.0 

SUMSQ-0.0 

ALPHA-XX 

N-N-l 
290  COF(2)-COF(2)+ALPHA*COF(l) 

DO  300  L-2,N 
300  COF ( L+l )-COF( L+l ) +ALPHA*COF( L)-SUMSQ*COF( L-l ) 
310  ROOTI(N2)-YY 

ROOTR(N2)-XX 

N2-N2+1 

IF  (SUMSQ)  320,330,320 
320  YY— YY 

SUMSQ-0.0 

GO  TO  310 
330  IF  (N)  30,30,90 

END 
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